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HIGHLIGHTS 


• First paper dealing with numerical dynamic simulation of thermal-lag Stirling engine. 

• Observation on different operating modes, such as rotating, swinging, swinging-to-rotate, and swinging-to-decay modes. 

• Study of influential factors affecting the operating modes. 

• Optimal geometrical or operating parameters of the baseline case for maximum engine power are found. 
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The present study is concerned with dynamic simulation of thermal-lag Stirling engines. A dynamic 
model is built and incorporated with a thermodynamic model to study the engine start process. A pro¬ 
totype engine is designed and simulated by using the dynamic model. In the simulation, different oper¬ 
ating modes, including rotating mode, swinging mode, swinging-to-rotate mode, and swinging-to-decay 
mode, have been observed. The rotating mode is desired and can be achieved if the operating parameters 
are properly designed. In a poor design, the engine may switch to the swinging or even the swinging-to- 
decay mode. In addition, it is found that geometric parameters, such as bore size, stroke, and volume of 
working spaces, also determine the operating mode of the engine. Brake thermal efficiency of the engine 
is monotonically reduced by increasing engine speed. However, study of the dependence of the shaft 
power of the engine speed shows that there exists a maximum value of the shaft power at an optimal 
operating engine speed. The optimal engine speed leading to maximum shaft power is significantly influ¬ 
enced by the geometrical parameters. 

© 2013 Elsevier Ltd. All rights reserved. 


1. Introduction 

Among the promising renewable energy technologies, dish con¬ 
centrating solar power (DCSP) systems have received increasing 
attention from the researchers in recent years. The DCSP system 
uses a dish to concentrate a large area of sunlight or solar thermal 
energy onto a small area. Electrical power is produced when the 
concentrated sunlight is converted to high-temperature thermal 
energy to drive a Stirling engine that is connected to an electrical 
power generator [1-3]. Besides, Stirling engine can also use natural 
gas or gasoline as auxiliary heat source so that it may still delivery 
power whenever the sunlight is not available. In some other cases, 
the engine is applied as a micro combined heat and power gener¬ 
ation (micro-CHP) unit [4,5] that can supply both residential hot 
water and electricity. 
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Earlier technological development in the Stirling engines was 
reviewed by Thombare and Verma [6]. Traditionally, Stirling en¬ 
gines are categorized into three types namely, oe-, p- and y-types, 
in terms of the arrangement of pistons, displacers, and cylinders. 
In the three types of engines, typically an engine is equipped with 
at least two moving parts, a piston plus a displacer or another pis¬ 
ton, in the cylinders. Relative movement between the two moving 
parts inside the cylinders is so carefully designed by assigning a 
phase angle between the two moving parts that the working fluid 
is able to flow back-and-forth between the heating and the cooling 
parts of the engine to absorb or reject heat, respectively, and then 
complete a thermodynamic cycle. Cheng and Yang [7] discussed 
the major differences in the configurations among the three types 
of engines and investigated relative performance of them. Rochelle 
and Grosu [8] took into consideration effects of exo-irreversibility 
and imperfect regeneration in the Schmidt model and obtained 
analytical solutions for engine speed, power and efficiency for dif¬ 
ferent types of Stirling engine. 

In 1979, Ceperley [9] presented a pistonless thermoacoustic 
Stirling engine which was demonstrated by using a looped tube 
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Nomenclature 

a A , a P acceleration of joints A and P (m/s 2 ) 
a Py y component of a P (m/s 2 ) 

a PA C gx. apAcg y x and y components of acceleration of center of 
_ gravity of link PA (m/s 2 ) 

AP vector (m) 

A r cross section void area of regenerative heater (m 2 ) 

c p constant pressure specific heat (J/kg K) 

c v constant volume specific heat (J/kg K) 

d p diameter of working space and piston (bore size) (m) 

D h hydraulic diameter of regenerative channel (m) 

e x , e y components of vector AP (m) 

/ Fanning friction factor 

fgt, f tg direction factor 

Fax. F Ay , F Px , Fpy x and y components of forces acting on joints A 
and P(N) 

Ff friction force at interface between piston and cylinder 

(N) 

G acceleration of gravity (m/s 2 ) 

h r heat transfer coefficient of regenerative heater (W/ 

m 2 K) 

If, / PA moment of inertia of flywheel and link PA (kg m 2 ) 

i. j. k unit vectors in Cartesian coordinate 
k conductivity of working fluid (W/m I<) 

Kf empirical coefficient in relation between velocity and 

pressure drop 

/ PA , l r length of link PA and regenerative heater (m) 
m g1 m t mass of working fluid in gas chamber and thermal buf¬ 
fer chamber (kg) 

nip, nip a mass of piston and link PA (kg) 

Nu Nusselt number in regenerative heater 

OA vector (m) 

p b , p g , p t pressure in buffer, gas, and thermal buffer chambers 
(Pa) 

Po. Pbo charged pressure in working space and buffer chamber 
(Pa) 

P perimeter of regenerative heater channel (m) 

Q in total heat input rate (W) 

Qi„g, Qj„ r, Qi n ,t rate °f heat transfer in gas chamber, regenerative 
heater, and thermal buffer chamber (W) 
r offset distance from crank to shaft center (m) 

R gas constant (J/I< mol) 

Re Reynolds number of flow through regenerative heater 


R c , R g , Rt thermal resistance in compression chamber, gas cham¬ 
ber, and thermal buffer chamber (K/W) 

Ri total thermal resistance in cold side (K/W) 

s coefficient in Eq. (18) (m -1 ) 

t time (s) 

T temperature of working fluid through regenerative hea¬ 

ter (K) 

T g , T t temperature of working fluid in gas chamber and ther¬ 
mal buffer chamber (K) 

T' g ,T' t temperature of working fluid at outlet of regenerative 
heater in gas chamber and thermal buffer chamber (K) 
T wg , T wt wall temperatures of gas chamber and thermal buffer 
chamber (K) 

Fwr temperature of regenerative porous matrix (K) 
u average velocity of working fluid through regenerative 

heater (m/s) 

v A , v P velocity of joints A and P (m/s) 
v Py y components of v P (m/s) 

V c , V g , V t volume of compression chamber, gas chamber, and 
thermal buffer chamber (m 3 ) 

VVj, VV S indicated power and shaft power (W) 

Greek symbols 

y c p /c v 

r f , r o torque due to friction force and working gas pressure 

(N m) 

r, oai loading torque (N m) 

>lt, >Im, >lb thermal, mechanical, and brake thermal efficiencies 
p dynamic viscosity of working fluid (N s/m 2 ) 

p p coefficient of kinetic friction 

(/> crank angle (rad) 

peg density of working fluid through regenerative heater 
(kg/m 3 ) 

t period (s) 

o> 0 , copa angular velocity of flywheel and link PA (rad/s) 

m a ve average angular velocity of flywheel (rad/s) 

co 0 , cijp A angular acceleration of flywheel and link PA (rad/s 2 ) 

Superscripts 
0 initial 

k time step 


with a differentially heated regenerator. The prototype engine with 
no moving piston revealed a concept that the traveling-wave ther¬ 
moacoustic oscillations can be regarded as a kind of heat engine. 
Lately, Swift [ 10,11 ] and Swift and Garrett [12] found that the ther¬ 
moacoustic concepts could be equally applicable to some recipro¬ 
cating heat engines. Hamaguchi et al. [13] developed a simple 
engine called “pulse tube engine”, which is equipped with one sin¬ 
gle reciprocating piston and does not utilize the natural frequency 
as in the thermoacoustic Stirling engines. Yoshida et al. [14] mea¬ 
sured work flux density distribution over the cross section of the 
pulse tube to elucidate the work source of the engine. They found 
that the pulse tube engine belongs to the standing wave engine 
group and the work source resides not in the porous stack but in 
the pulse tube. 

On the other hand, one other concept of the Stirling engine was 
proposed by Chen and West [15] and the engine was entitled 
“thermal-lag Stirling engine” afterward. The thermal-lag Stirling 
engine requires no displacer and its connected link, and instead, 
a porous-medium stack is fixed in the cylinder as a static regener¬ 
ative heater. The first prototype of the thermal-lag engine and its 


performance testing were made by Tailer [16,17]. He stated that 
the engine experiences imperfect cooling and heating processes 
in the working spaces, and the imperfect cooling and heating cause 
a time lag in the gas temperature response as following the move¬ 
ment of the piston. Thus, the engine was named “thermal-lag” en¬ 
gine. A numerical model for the thermal-lag engine was proposed 
by Arques [18]. In this report, it could be observed that the areas of 
the thermodynamic cycles of the thermal-lag engine plotted in the 
p-V diagrams are rather small. It seemingly implies that that the 
thermal-lag engines are not capable of producing appreciable work 
output. Recently, Cheng and Yang [19] made a prototype thermal- 
lag engine of 15 W and also developed a thermodynamic model for 
predicting transient variations of the thermodynamic properties in 
the individual working spaces of the engine. Effects of geometrical 
and operating parameters, such as heating and cooling tempera¬ 
tures, volumes of the chambers, thermal resistances, stroke of pis¬ 
ton, and bore size on indicated power output and thermal 
efficiency were also evaluated. 

Similar to the pulse-tube engine, the thermal-lag engine has 
only one moving part (piston) and a static part (regenerative 
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heater) without the displacer and its link. As a matter of fact, the 
thermal-lag engine is frequently referred to as a special kind of 
the pulse tube engines. However, a standard pulse tube engine 
consists of a piston in one cylinder, a pulse tube, an orifice, a hea¬ 
ter, a regenerator, a cooler and a second cooler, and a buffer zone 
[10-14], whereas the mechanism of the prototype of the ther¬ 
mal-lag engine presented by Tailer [16,17] does not include the 
pulse tube, the orifice, the second cooler, and the buffer zone so 
that its mechanical structure can be further simplified. Due to its 
simplicity in structure, its cost in manufacturing and maintenance 
can relatively low. Besides, it has been observed from the experi¬ 
ments on the prototype engine made by Cheng and Yang [19] that 
the engine can be slightly pushed to start. When a small initial 
rotation speed is applied, the engine starts with a bi-directional 
swinging mode in a short period of time and then proceeds to a 
uni-directional rotating mode. This can be called “swinging-to- 
rotate” mode. It implies that the engine is possible to start auto¬ 
matically if engine friction is improved. These advantages make 
this engine competitive particularly when it is applied in a small 
scale DCSP system for converting solar energy into electricity in 
the mountains or other remote areas where regular maintenance 
is difficult. Furthermore, in addition to the swinging-to-rotate 
mode, authors observed other different operating modes of the 
prototype engine in the experiments, and the physical and geomet¬ 
rical factors affecting the operating modes are still not clear. Unfor¬ 
tunately, no numerical model is available for the dynamic 
simulation of the start process of the thermal-lag engine. There¬ 
fore, in this study a dynamic model is built and incorporated with 
the thermodynamic model to study the engine start process and 
clarify the influential factors affecting the operating factors. 

The schematic of the thermal-lag Stirling engine is plotted in 
Fig. 1. The space of the engine is separated into gas chamber, 
regenerative heater, thermal buffer chamber and compression 
chamber by the regenerative heater and the piston. As described 
earlier, the piston is the only one moving part reciprocating be¬ 
tween top-dead and bottom-dead points. Therefore, the volumes 
of gas chamber, regenerative heater chamber and thermal buffer 
chamber are fixed, whereas the volume of the compression cham¬ 
ber is varied with time. Wall temperatures of the gas chamber and 
compression chamber (hot part) and the thermal buffer chamber 
(cold part) are maintained at T wg and T w[ , respectively. An axial 
temperature gradient is hence established within the regenerative 
heater. The regenerative heater is made of metal wires, screen 
mesh, or metal foams, which serves not only for heating but also 
for resulting in a phase angles in the periodic variations of volumes 


Table 1 

Geometrical and operating parameters of base-line case. 


Geometrical parameters 

r (m) 

0.025 

/pa (m) 

0.075 

dp(m) 

0.045 

V, (m * 1 2 3 4 ) 

5 x 10" 6 

Mm 3 ) 

30 x 10"‘ 

m p (kg) 

0.055 

m PA (kg) 

0.055 

7pa (kg m 3 ) 

5 x 10“ 5 

7/(kg m 3 ) 

0.005 

Operating parameters 

cu° (rpm) 

500 

T w , (K) 

300 

T wg (K) 

1000 

R g (K/W) 

1 

R, (K/W) 

200 

Rc (K/W) 

1 

Po = Pbo (Pa) 

101,325 

Pp 

0.8 


and pressures. The prototype engine is regarded as the baseline 
case in this study. The geometrical and operating parameters of 
the prototype engine are listed in Table 1. According to experimen¬ 
tal observation on the engine, there exist four possible operating 
modes in the start process: 

(1) Rotating mode: the flywheel starts to rotate immediately 
and the rotation speed smoothly increases to steady-state 
regime. 

(2) Swinging mode: the flywheel exhibits a stable swinging 
motion but without rotation. 

(3) Swinging-to-rotate mode: the flywheel starts with a swing¬ 
ing mode in a short period of time and then proceeds to a 
uni-directional rotating mode. 

(4) Swinging-to-decay mode: the flywheel starts with a swing¬ 
ing mode then the swinging motion gradually decay to a full 
stop. 

The swinging motion is a unique behavior with the thermal-lag 
Stirling engines, which was also observed by Tailer [16,17]. It is 
recognized that the operating mode is dependent not only on the 
geometrical parameters but also on the operating conditions. 

2. Dynamic model 


Gas chamber 




Fig. 1. Schematic of thermal-lag Stirling engine with crank drive mechanism. 


In this section, the dynamic model for predicting the dynamic 
behavior of the thermal-lag Stirling engine is described briefly. 
The photo of the prototype engine is shown in Fig. 1. The geometric 
parameters of the thermal-lag Stirling engine with crank drive 
mechanism are also illustrated in this figure. The crank drive 
mechanism consists of three parts, namely piston, link, and crank, 
which is connected to the flywheel. Free body diagrams of the parts 
are shown in Fig. 2. Pressure force of the working gas is transferred 
from the piston to the flywheel to push the flywheel rotating. The 
forces exerted on the parts at any time step can be determined 
based on force and moment balance principles. Therefore, crank 
angle, angular speed and acceleration of the flywheel are deter¬ 
mined with the help of the dynamic model, and positions, veloci¬ 
ties and accelerations of the parts can be readily calculated based 
on the kinematics analysis [20]. According to the information of 
positions, velocities and accelerations of the parts, one can deter¬ 
mine the volumes of all the chambers in the engine, and then pre¬ 
dict the pressures in the chambers for the consecutive time step 
with the help of the thermodynamic model. 

In the dynamic model, the following assumptions are made: 
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(b) Piston 



(c ) Link PA 



Fig. 2. Free body diagrams of crank drive mechanism. 


1. All moving parts are rigid bodies with no deformation. 

2. Magnitude of friction force between the piston and the cylinder 
is proportional to the normal force to the contact interface. 

3. Friction force at each joint is negligible. 

4. Center of gravity of each moving part is at its geometric center. 
2.3. Velocities and accelerations 

In accordance with Fig. 2, the velocity and acceleration of joint A 
are: 

v A = fflo x 0A k = -C 0 or(sin<j> k i - cos <t> k j) (1) 

a A = chg x 0A k + ffig x x OA k j 
= - ^cbgrsin <j> k + (a>o) 2 rcos ^ k J i 
+ |d>orcoS(/) k - (cOo) 2 rsin<j> k J j (2) 


where <j>\ = ccfc k, and cb^ = to £ k are crank angle, angular veloc¬ 

ity, and angular acceleration, respectively, at time step k in itera¬ 
tion. And, the velocity of piston can be calculated by: 

vii = vi + mjS, x AP k 

_(3) 

= sin <j> k - <u|S A e k ) i + (co^r cos <j> k - to |S A e k ^) j 

where e k = -rcos<j> k and e k = \Jlp A - r 1 2 3 4 cos 2 <j> k representing the x- 

and y-components of vector AP. Note that x-component of v k is 
equal to zero. One obtains angular velocity of link PA and y-compo- 
nent of v k as 


ojpA = -oo^r sin <j> k /e k 

( 4 ) 

and 


v k y = aJo r ( cos <t> k - e x sin <t> k / e k ) 

( 5 ) 


The acceleration of joint P is determined by 

ap = a A + cbp A x AP k + ft)p A x ( a)p A x AP k ) 


( 6 ) 
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Note that x-component of a£ is equal to zero and and y-compo- 
nent of a£ can be expressed as 

“pa = - [“o rsin / + (oJq ) 2 r cos <j> k + (<u k A ) 2 e k ]/e k (7) 

and 

a k y = [cb^r cos <4 - (m£) 2 r sin <4 - (<4 A ) 2 e k ] 

- [cb k rsin<4 + (<u k ) 2 rcos <4 + (co k m ) 2 e k ^e k /e k y (8) 


(1) Working fluid is air, which is treated as an ideal gas and con¬ 
tained in the engine without leakage. 

(2) Pressure inside the engine is uniform. 

(3) Wall temperature of the gas chamber is maintained at T wg 
due to heating by an external high-temperature reservoir, 
and on the other hand, that of the thermal buffer chamber 
and the compression chamber is maintained at T w t due to 
external cooling. 

(4) Average flow velocity in the regenerative heater is propor¬ 
tional to the pressure drop. 


2.2. Equations of motion 

Let F k x and F k y denote the x- and y-components of the force act¬ 
ing on the piston by link PA, and p k and pj[ the pressures of the 
working gas exerted on the top and bottom surfaces of the piston, 
respectively. One has the equation of motion for the piston as 

-Tid 2 (p k - p\)/A - m p g + Fp y = m„a k y (9) 

F k y can be determined by the above equation. 

According to force and angular momentum balance, the 
momentum and the angular momentum equations for link PA 
can be written as: 

m PA a PAcgx = - f Px ( 10a ) 

mpAClpAcgy = “ F Py ~ m PAg (10b) 


Herein the gas chamber is called hot side, and the combined 
space of the thermal buffer chamber and the compression chamber 
is called cold side. An axial temperature gradient is hence estab¬ 
lished within the regenerative heater. Mass flow rate through the 
regenerative heater is determined in terms of the pressure differ¬ 
ence between the hot and the cold sides and the properties of 
the porous medium as well. For a regenerative heater which is 
composed of metal form or wire-mesh woven-screen matrices, 
the pressure drop can be expressed in dimensionless form as a Fan¬ 
ning friction factor / [21], which is defined as / = where 

lr P tg U 

A p = pi - p k and hydraulic diameter D h = 4 A r /P. 

The form of the relation between the Fanning friction factor and 
the Reynolds number corresponding to fully developed laminar 
flows in micro channels [22] is adopted here as: 

/Re = Kf (15) 


/pa«4 = (4$ - F k P X + 4*4 - F *y4) /2 (11) 

One can solve Eqs. (10) and (11) forF k x , F^, and F Ay and then obtain 
the torque on joint 0 as 

44 = r ( F L sin 4> k ~ 4y cos (,b k ) (12) 


The angular acceleration of the flywheel ebo +1 caused by 4 +1 , rf +1 
and a loading torque is determined by the equation of angular 
motion for the flywheel: 


0 k+l 


( pk+1 I r’k+l 

{* 0 + 1 f 



(13) 


where If is moment of inertia of the flywheel, and Jy+J is deter- 
mined by the operating conditions, and F k+1 torque loss due to fric¬ 
tion which is calculated by r k+1 = F k r cos <j> k . F k is friction force at 
the interface between the piston and the cylinder determined by 

4 = ~/'p| 4*|4/|4|- 


Integrating Eq. (13) with respect to time, one yields the angular 
velocity and the crank angle at the consecutive time step. As a re¬ 
sult, the volumes of the compression and the buffer chambers can 
be obtained and introduced into the thermodynamic model, which 
is illustrated in the subsequent section. 

The obtained results of the instantaneous crank angle and 
angular velocity co 0 of the flywheel can be used to estimate the 
average angular velocity in a cycle as: 


tttave 


1 

T 



(14) 


3. Thermodynamic model 


where Kf is an empirical coefficient. Relation between average 
velocity and the pressure drop can then be derived as 

u k = -D 2 h (p k g -pl)/2pK f l r (16) 

The mass flow rate through the regenerative heater can be calcu¬ 
lated by m k = pl g A r u k . The average velocity is assigned to be posi¬ 
tive when the working fluid flows from the thermal buffer 
chamber to the gas chamber. 

The axial temperature profile of the regenerative porous matrix is 
assumed to be a steady and linear function varied from the hot to the 
cold sides. When the working fluid passes through the regenerative 
heater, the enthalpy change of the working fluid is mainly caused by 
the convective heat transfer to or from the working fluid. Therefore, 

pl g c p u k A r dT = h r Pdx(T wr - T) (17) 

where T wr is the temperature of the regenerative porous matrix, and P 
is the cross section void perimeter of the regenerative porous matrix. 
Since the axial temperature profile of the regenerative porous matrix 
is assumed to be a steady and linear function varied from T wg to T w[ , 
one has T, = — Twt j Tws x + T wt . As the working gas flows from the cold 
to the hot sides, the working gas absorbs heat from the regenerative 
heater and hence its temperature rises from T k to T*. On the contrary, 
while the working gas flows from the hot to the cold sides, it is cooled 
and its temperature falls from T k to T g . Therefore, the boundary con¬ 
ditions associated with Eq. (17) is described as: 

(1) If u k > 0, T = T k at x = 0: 

One obtains the temperature distribution of the working fluid in 
the regenerative heater as: 

T(x) = (T k - T wt )e- k * - + T wt (18) 


The thermodynamic model presented by Cheng and Yang [19] is 
extended and is incorporated with the dynamic model. To save the 
space of the report, the thermodynamic model is described briefly 
in this section. 

In thermodynamic model, some assumptions were made as 
follows: 


where s k = 4kNu/p k g c p u k D^, and k is thermal conductivity of the 
working fluid. Let x = l r , and then Eq. (18) gives the temperature 
of working fluid at the exit of the regenerative heater to the hot side 

T g = T wg - jq^(F wg - T wc ) + (T k - T wt )e^ (19a) 
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(2) If u k < 0, T = T k at x = 0: 

Similarly, one may obtain the temperature of working fluid at 
the exit of the regenerative heater to the cold side 

m _jk| 

r t k = T wt + (T wg -T wl )~ ( T wg - T k ) e-* 1 ' (19b) 

Transient temperature variation of the working fluid in the cold side 
(thermal buffer and compression chambers) can be predicted based 
on the energy conservation law. In finite-difference form, it is ex¬ 
pressed as 


yk+1 

* t 


T k Pt A/k+I 
m k+l t m M Cv V c 


+ 


At 

m't+'c,, 



( 20 ) 


On the right-hand side of the above equation, the second term is 
the work done by the working fluid, and the first term in bracket 
denotes the heat transfer which is calculated in terms of thermal 
resistance and temperature difference. In addition, the second 
term in the bracket represents energy transfer rate accompanying 
the mass flow, and hence, it is dependent on the direction of the 
mass flow. The two direction factors f k g and j J appearing in Eq. 
(20) are introduced to show the direction effect, which are defined 
as 


k = J 1 if u k > 0 k 

Jtg \0 if u k < 0' J * 


1 if u k <0 
0 if u k > 0 


( 21 ) 


Similarly, temperature of the working fluid in the cold side (gas 
chamber) can be calculated by: 


yk+1 

g 


m s T k A t 

m k+i s mjt+’Ct, 


r "s T g ! ^k f 

pkf/k J_ fk-rk) 

R g +mC ?( J 

tg l g ^Jgt^g) 


( 22 ) 


In Eqs. (20) and (22), R ( and R g represent the thermal resistances in 
the cold and the hot sides, respectively. R/ is calculated by [19] 


R 


k 

1 


1 

R t 


l 

R c 




(23) 


There is no moving-boundary work done by the working fluid in the 
hot side because the volume of the gas chamber is fixed. With the 
ideal-gas equation of state, the pressures in the cold and the hot 
sides are calculated, respectively, as 


P 


k+l 


m k+1 RT k+ 1 

v k+1 + V t 


(24a) 


P 


k+l 

g 


m^KT 


k+l 

g 


Vg 


(24b) 


Total heat input and indicated power can be calculated by: 

Q.in=\f^ M + Q.%n + Q.lindt (25a) 


and 


W, = 1 j p k dV k c (25b) 

where Q k in , Q k in , and Q k in are the heat transfer rates in the thermal 
buffer chamber, the regenerative heater, and the gas chamber, 
respectively, and p k is the pressure in the thermal buffer and the 
compression chambers. In addition, the thermal efficiency (»/ t ) can 
be obtained by 


+ = Wi/Q.in (26) 

As a loading torque I ,k oad is applied on the engine shaft, the output 
shaft power of the engine can be yielded by 

Ws = \fr k oad co k dt (27) 

Thus, the mechanical efficiency (f/ m ) and the brake thermal effi¬ 
ciency (r]t) can be obtained in terms of the works based on the fol¬ 
lowing definitions: 

1m = W,/W, (28) 

and 



Fig. 3. Flow chart of computation process. 
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’lb = ’h’lm = W s /Q_ in (29) 


The pressure, the temperature, the instantaneous heat transfer rate 
and the working gas mass in each of the chambers can be deter¬ 
mined with the help of the thermodynamic model. The information 
of the pressures in the chambers is then introduced to the dynamic 
model to calculate the crank angle, the angular velocity, and the 
angular acceleration of the flywheel iteratively. The flow chart of 
computation process is shown in Fig. 3. As has been described 
above, the dynamic behavior of thermal-lag Stirling engine can be 
simulated by combining the dynamic and the thermodynamic mod¬ 
els. The geometrical parameters and operating conditions of the 
baseline case are listed in Table 1. The initial crank angle and angu¬ 
lar acceleration are both assigned to be zero. The initial pressures in 
all the chambers are set to be equal to the atmospheric pressure. 
The temperatures of the hot and cold parts are maintained at 300 
and 1000 K, respectively. 



Fig. 4. Instantaneous and average angular velocities. 



Fig. 5. Instantaneous angular velocity under different initial angular velocities. 


4. Results and discussion 

For the base-line case, the instantaneous and the cyclic average 
angular velocities, co 0 and to ave , are plotted in Fig. 4. The engine 
speed increases from the initial angular velocity of 500 rpm to 
steady state regime at 1980 rpm in 100 s, and the fluctuation of 
instantaneous angular velocity is about 75 rpm. Compared with 
traditional Stirling engine, the present thermal-lag Stirling engine 
needs longer time (say, 10 times that with the traditional engines) 
to reach the steady-state regime. This is because there is no dis¬ 
placer in the thermal-lag Stirling engine to move the working gas 
immediately. In this regard, it needs more time to adjust the phase 
between the pressure and the volume variations in working spaces. 

Fig. 5 shows the instantaneous angular velocities under differ¬ 
ent initial angular velocities: (Mq = 500 rpm, ct>o = 1500 rpm, and 
cOq = 2500 rpm. It is expected that instantaneous and the cyclic 
average angular velocities eventually lead to equal cyclic average 
angular velocity in the steady-state regime. In this figure, all the 
curves of start processes associated with different initial rotation 
speeds exhibit smooth changes in the rotation speed from start 
to steady-state regime. This is called rotating mode in this study. 



(a) Swinging-to-rotate mode 



Fig. 6. Effects of low initial angular velocity. 
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Time (sec) 

(a) Rotating mode 



Time (sec) 


(b) Swinging mode 

Fig. 7. Effects of moment of inertia of flywheel. 

Fig. 6a shows the instantaneous angular velocities at two initial 
rotation speeds, 150 and 160 rpm. As the initial angular velocity is 
lowered, the flywheel swings first till the amplitude of the piston 
oscillation becomes lager than a half stroke, and then the flywheel 
starts to rotate randomly in either clockwise or counter-clockwise 
direction. This is referred to as the swinging-to-rotate mode. It is 
also found that no matter which direction the flywheel chooses, 
eventually the magnitude of engine speed will be the same. At a 
lower initial angular velocity, the flywheel spends longer time in 
the swinging stage before the steady-state regime is achieved. 
When the initial angular velocity is further lowered to 145 rpm, 
the flywheel still starts with a swinging mode; however, the 
swinging motion gradually decays to a full stop, as illustrated in 
Fig. 6b. Here one observes the swinging-to-decay mode. 

Fig. 7 shows the effects of the moment of inertia of the flywheel 
on the transient variation in the instantaneous angular velocity. In 
Fig. 7a, the cases at //= 0.005, 0.03, and 0.05 kg m 2 are considered, 
and the initial angular velocity is fixed at ct>^ = 500 rpm. It is seen 
that the operating mode with all the cases shown in Fig. 7a are 
rotating mode. With a larger moment of inertia, it takes longer to 
reach the steady-state regime, and nevertheless the fluctuation of 



Time (sec) 

Fig. 8. Instantaneous angular velocity variation under different loading torques. 

the instantaneous angular velocity is smaller. The average angular 
velocity in the steady-state regime is independent of the moment 
of inertia. 

When the moment of inertia of flywheel is lowered to 
0.0001 kg m 2 , it is interesting to find that the angular velocity 
oscillates stably between positive and negative values, as depicted 
in Fig. 7b. This actually implies that the flywheel is swinging but 
not able to rotate. That is called the swinging mode in this study. 
In the swinging mode, the amplitude of piston oscillation is always 
less than a half stroke, and thus the engine is unable to turn out to 
the rotating mode. According to the obtained results, it can be con¬ 
cluded that the operating and geometrical parameters indeed have 
profound influence on the operating mode of the engine. 

Fig. 8 shows the variation in instantaneous angular velocity for 
the base-line case under different loading torques. The engine is 
started without loading to the steady-state regime in the first 
100 s, and then the loading torque (-T| 0a d) is applied. As a result, 
the average angular velocity decreases remarkably and attains an¬ 
other steady-state regime. The average angular velocity in the lat¬ 
ter steady-state regime is reduced if the loading torque is elevated. 
It is also observed that the fluctuation of instantaneous angular 
velocity becomes larger with higher loading torque. 

Fig. 9 shows the effects of heating temperature and thermal 
resistance on shaft power and brake thermal efficiency. In Fig. 9a, 
the effects of the heating temperature on the performance of the 
engine at various heating temperatures are evaluated. It is ob¬ 
served the both the shaft power and brake thermal efficiency in¬ 
crease with the heating temperature. At a heating temperature 
within 800-1200 K, the shaft power reaches a maximum at about 
500 rpm. It indicates that the optimal rotation speed of the proto¬ 
type engine is 500 rpm to gain a maximum power. It is important 
to mention that the thermal efficiency of the thermal-lag engines 
are relatively low compared to other types of Stirling engines be¬ 
cause with no displacer the working fluid might not be swept 
effectively between the hot and cold parts. In Fig. 9a, for all cases 
the brake thermal efficiency is below 14%. Therefore, there is still 
room to increase the performance of the thermal-lag engines. 
Fig. 9b conveys the performance of the engine with various ther¬ 
mal resistances. It is clearly observed that to increase the shaft 
power and the brake thermal efficiency of the engine, the thermal 
resistances in the hot and the cold parts ( R g and R c ) should be re¬ 
duced. For the case at R g = R c = 0.5 K/W, the maximum shaft power 
is about 30 W at 1000 rpm. However, at R g = R c = 1.0 K/W, the max- 
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CO ave (rpm) 

(b) Effect of thermal resistance 


Fig. 9. Effects of heating temperature and thermal resistance on shaft power and 
brake thermal efficiency. 


imum shaft power is decreased to 15 W. It is important to mention 
that if the thermal resistance is too large, the operating mode of the 
engine may switch from the rotating mode to the undesired swing- 
ing-to-decay mode. Fig. 10 shows the results for the case at R g = - 
R c = 10K/W. In this figure, all other parameters are fixed at the 
same values as the baseline case except the thermal resistances. 
The rotation speed is varied from start to the swinging mode at 
such a higher thermal resistance. This may be attributed to insuf¬ 
ficient heat transfers in the hot and cold parts caused by higher 
thermal resistances. 

Effects of stroke and bore size on the shaft power and the brake 
thermal efficiency are displayed in Fig. 11. The numerical data 
shown in Fig. 11a indicate that the case with smaller stroke has 
higher average angular velocity. Nevertheless, with a larger stroke, 
the compression ratio is higher and heat rejection is more suffi¬ 
cient; therefore, one may have higher shaft power and brake ther¬ 
mal efficiency. Meanwhile, a larger stroke causes larger friction 
loss. Therefore, there is an optimal stroke for maximum shaft 
power. In this figure, the optimal stroke is r = 0.03 m and the cor¬ 
responding maximum shaft power is 16.95 W at 415 rpm. In 



Fig. 10. Rotation speed (from start to swinging mode) at R g = R c = 10 1</W. 



(a) Effect of stroke 



(b) Effect of bore size 

Fig. 11. Effects of stroke and bore size on shaft power and brake thermal efficiency. 
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Fig. 12. Rotation speed (from start to swinging mode) at d p = 0.07 m. 



(a) Effect of thermal buffer chamber 



(b) Effect of gas chamber 


Fig. 13. Effects of volumes of thermal buffer chamber and gas chamber on shaft 
power and brake thermal efficiency. 


Fig. lib, the effects of bore size are displayed. It is found that a case 
with larger bore size has lower average angular velocity and needs 
a higher initial rotation speed to start the engine. However, on the 
other hand, with a larger bore size, the cylinder contains more 
mass of the working gas; therefore, the engine may produce higher 
shaft power. Hence, one may find an optimal bore size for maxi¬ 
mum shaft power. In this figure, it is observed that the optimal 
bore size is d p = 0.05 m, and the corresponding maximum shaft 
power is 17.69 W at 411 rpm. 

If the bore size is too large, the engine may switch to the swing¬ 
ing mode or even swinging-to-decay mode. As plotted in Fig. 12, at 
d p = 0.07 m, the engine starts from initial instantaneous angular 
velocity of 500 rpm and eventually attains the swinging mode. In 
this situation, the flywheel is not able to rotate. 

Fig. 13 shows the effects of the volumes of the thermal buffer 
chamber and gas chamber on the performance of the engine. In 
Fig. 13a, it is found that a reduction in the volume of the thermal 
buffer chamber leads to higher performance. This is because the 
thermal buffer chamber is regarded as a dead space in the ther¬ 
mal-lag Stirling engine which should be minimized. Fig. 13b de¬ 
picts the performance curves of the engine for various volumes 
of the gas chamber. It is seen that there exists an optimal gas 
chamber volume for a maximum shaft power. In this figure, the 
optimal volume of the gas chamber is found to be 1^=10 cm 3 , 
and the corresponding maximum shaft power is 17.8 W at 
587 rpm. However, if a maximum brake thermal efficiency is of 
interests, the optimal volume of the gas chamber would be 
V g = 5 cm 3 . 

5. Conclusions 

Numerical dynamic simulation of the thermal-lag Stirling en¬ 
gines by incorporating a dynamic model with a thermodynamic 
model has been performed. In the simulation, different operating 
modes, including rotating mode, swinging mode, swinging-to-ro- 
tate mode, and swinging-to-decay mode, have been observed. 
The rotating mode is desired and can be achieved if the operating 
and geometrical parameters are properly selected. In a poor design, 
the engine may switch to the swinging or even the swinging-to-de- 
cay mode. For the particular cases considered here, the swinging 
mode is observed for the cases with a larger bore size 
(d p = 0.07 m) or larger thermal resistances (R g = R c = 10 K/W). 

It is found that the brake thermal efficiency of the engine is 
monotonically reduced by increasing engine speed. However, para¬ 
metric study of the dependence of the shaft power of the engine 
speed shows that there exists a maximum value of the shaft power 
at an optimal operating engine speed. There also exist optimal 
stroke, bore size, and volume of gas chamber for achieving maxi¬ 
mum shaft power or higher brake thermal efficiency. Through 
the parametric study for the baseline case, it has been found that 
the optimal stroke is r = 0.03 m, optimal bore size d p = 0.05 m, 
and optimal volume of the gas chamber V g = 10 cm 3 . 
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